Install and load ciftiTools:
if (!require("ciftiTools", quietly=TRUE)){
install.packages("ciftiTools")
} else if (packageVersion("ciftiTools") < "0.9.0") {
update.packages("ciftiTools")
}
library(ciftiTools)
Download the Connectome Workbench here: https://www.humanconnectome.org/software/get-connectome-workbench
Tell ciftiTools where to find the Connectome Workbench:
# Replace '~/Desktop/workbench' with the actual path on your computer.
# If successful, the path to the Workbench executable will be printed.
ciftiTools.setOption("wb_path", "~/Desktop/workbench")
## Using this Workbench path: '/Users/ddpham/Desktop/workbench/bin_macosx64/wb_command'.
Setup rgl plots for knitting.
library(rgl)
## Warning: package 'rgl' was built under R version 4.1.2
rgl::setupKnitr()
# Sometimes the first OpenGL window does not render properly.
rgl::rgl.open(); rgl::rgl.close()
Read in our first CIFTI file!
cii_fname <- file.path("data", "MyelinAndCorrThickness.dscalar.nii")
xii <- read_cifti(cii_fname)
xii # same as `summary(xii)`
## ====CIFTI METADATA===================
## Intent: 3006 (dscalar)
## - names "MyelinMap_BC_decurv", "corrThickness"
## Measurements: 2 columns
##
## ====BRAIN STRUCTURES=================
## - left cortex 30424 data vertices
## 2068 medial wall vertices (32492 total)
##
## - right cortex 30527 data vertices
## 1965 medial wall vertices (32492 total)
Plot it:
# same as view_xifti(xii) and view_xifti_surface(xii)
plot(xii, zlim=c(1,2))
Myelin
View a list of commonly-used functions:
if(interactive()) { help(ciftiTools) }
View details for a specific function:
if(interactive()) { help(view_xifti_surface) }
Use read_cifti to read in the CIFTI file named cii_fname. Set the brainstructures argument to include all brain structures.
cii_fname <- file.path("data/toy_fMRI.dtseries.nii")
# [Use `read_cifti` here]
xii <- read_cifti(cii_fname, brainstructures="all")
xii
## ====CIFTI METADATA===================
## Intent: 3002 (dtseries)
## - time step 1 (seconds)
## - time start 0
## Measurements: 500 columns
##
## ====BRAIN STRUCTURES=================
## - left cortex 3670 data vertices
## 332 medial wall vertices (4002 total)
##
## - right cortex 3674 data vertices
## 328 medial wall vertices (4002 total)
##
## - subcortex 31870 data voxels
## subcortical structures and number of voxels in each:
## Cortex-L (0), Cortex-R (0),
## Accumbens-L (135), Accumbens-R (140),
## Amygdala-L (315), Amygdala-R (332),
## Brain Stem (3472),
## Caudate-L (728), Caudate-R (755),
## Cerebellum-L (8709), Cerebellum-R (9144),
## Diencephalon-L (706), Diencephalon-R (712),
## Hippocampus-L (764), Hippocampus-R (795),
## Pallidum-L (297), Pallidum-R (260),
## Putamen-L (1060), Putamen-R (1010),
## Thalamus-L (1288), Thalamus-R (1248).
We will add very inflated surfaces prior to plotting.
xii <- add_surf(
xii,
load_surf("left", "very inflated"),
load_surf("right", "very inflated")
)
Use view_xifti_surface to plot the first timepoint.
Customize your plot by specifying some arguments! Here’s a few ideas:
titlezlim (min, max)colors to the name of a sequential palette from RColorBrewer. To see a list, use RColorBrewer::display.brewer.all().# [Use `view_xifti_surface` here]
view_xifti_surface(
xii, title="T=1",
zlim=c(5000, 10000), colors="PuRd"
)
First timepoint
Here’s a function to calculate the fALFF of a single time series vector.
#' Adapted from `ANTsR::alffmap`
falff <- function(x, TR=.72, f_low=.01, f_high=.1){
# Subtract linear trend
x <- lm(x~seq(length(x)))$residuals
# Compute Fourier decomposition.
y <- stats::spec.pgram(
stats::ts(x, frequency=1/TR),
taper=0, plot=FALSE
)
freq_mask <- (y$freq > f_low) & (y$freq < f_high)
# Compute proportion within frequency range
sum(y$spec[freq_mask]) / sum(y$spec)
}
Use apply_xifti to apply the falff function to each row of the data matrix in xii.
# [Use `apply_xifti` here.]
fxii <- apply_xifti(xii, 1, falff)
Plot the results using a color scale going from 0 (white) to 1 (blue), and save to a PNG file. For help on finding the relevant arguments, read help(view_xifti_surface).
# [Use `view_xifti_surface` here.]
view_xifti_surface(fxii, colors=c("white", "blue"), zlim=c(0, 1))
fALFF
# [Use `view_xifti_volume` here.]
view_xifti_volume(fxii, colors=c("white", "blue"), zlim=c(0, 1))
fALFF, subcortex
view_xifti_surface(fxii, colors="Blues", zlim=c(0, 1))
fALFF, alternate colors
Read the fMRI data and parcellation.
cii_fname <- file.path("data/MSC03_rest07.4k.dtseries.nii")
xii <- read_cifti(cii_fname)
dim(xii)
## [1] 7344 818
parc <- read_cifti("data/Kong2022_100Parc.dlabel.nii", resamp_res=4000)
dim(parc)
## [1] 8004 1
# rownames(parc$meta$cifti$labels[[1]])
table(c(as.matrix(parc))==0)
##
## FALSE TRUE
## 7219 785
Incompatible ROIs: we need to replace the medial wall mask with NA values.
xii <- move_from_mwall(xii)
Pick our seed and plot it.
# View(parc$meta$cifti$labels$parcels["Key"])
plot(parc, title="Parcellation")
Parcellation
the_seed <- 3 # sample(seq(100), 1)
parc_masked <- transform_xifti(parc, function(q){ifelse(q==the_seed, q, 0)})
plot(parc_masked, borders="black", title=paste0("Parcel #", the_seed))
Selected parcel
Now that we’ve identified the seed, let’s do this!
Use select_xifti to mask out high-motion timepoints. (Only keep the columns for which the corresponding entry of xii_tmask is TRUE.)
xii_tmask <- read.csv("data/MSC03_rest07_tmask.txt")
xii_tmask <- which(as.logical(xii_tmask[,1]))
# [Use `select_xifti` here.]
xii <- select_xifti(xii, xii_tmask)
dim(xii)
## [1] 8004 475
The below code chunk is already filled out. Here’s what it does:
seed_mask: Convert the parcellation "xifti" to a numeric vector of parcel indices for each vertex, and then convert that to a binary vector indicating whether the parcel is the seed parcel.parc_mean_ts: Convert the resting-state fMRI "xifti" to a matrix, subset the rows to select only vertices inside the parcel, and then take the mean across space.xii_sc: Compute correlation between each vertex time series, and the parcel mean time series.seed_mask <- c(as.matrix(parc)) == the_seed
parc_mean_ts <- colMeans(as.matrix(xii)[seed_mask,,drop=FALSE])
xii_sc <- cor(t(as.matrix(xii)), parc_mean_ts)
xii_sc <- newdata_xifti(xii, xii_sc)
Now we will convert the result back to a "xifti" with newdata_xifti, and then format is as a “dscalar” since the data no longer represents a time series.
xii_sc <- cor(t(as.matrix(xii)), parc_mean_ts)
# [Use `newdata_xifti` here.]
xii_sc <- newdata_xifti(xii, xii_sc)
# [Use `convert_xifti` here.]
xii_sc <- convert_xifti(xii_sc, "dscalar")
Write out the result with write_xifti.
seed_fname <- paste0("seed_cor_", the_seed, ".dscalar.nii")
# [Use `write_xifti` here.]
write_xifti(xii_sc, seed_fname)
## Writing left cortex.
## Writing right cortex.
## Creating CIFTI file from separated components.
Finally, plot the result!
view_xifti_surface(xii_sc, zlim=c(-.7, .7))
Seed correlation
Install the latest version of BayesfMRI from GitHub.
if (!requireNamespace("BayesfMRI", quietly=TRUE) || packageVersion("BayesfMRI") < "0.2.0") {
if (!requireNamespace("devtools", quietly=TRUE)) { install.packages("devtools") }
devtools::install_github("mandymejia/BayesfMRI", "2.0")
}
library(BayesfMRI)
Read in and format the event onsets file. Use the “left foot” task.
onsets <- paste0("data/MSC03_motor07_run0", seq(2), "_events.tsv")
onsets <- lapply(as.list(onsets), read.csv, sep="\t")
onsets <- lapply(onsets, function(q){list(as.matrix(q[q[,3]=="LFoot",seq(2)]))})
Run classical GLM for the right and left cortex separately.
bglmR <- BayesGLM_cifti(
paste0("data/MSC03_motor07_run0", seq(2), ".4k.dtseries.nii"),
brainstructures="right", surfR = load_surf("right", resamp_res=4000),
onsets=onsets, ar_order=0, EM=FALSE, Bayes=FALSE, resamp_res=NULL, TR=2.2
)
##
## SETTING UP DATA
## .. reading in data for session 1
## .. reading in data for session 2
## MAKING DESIGN MATRICES
## RUNNING MODELS
##
## .. RIGHT CORTEX ANALYSIS
## 328 locations removed due to NA or NaN values in at least one scan.
## PUTTING RESULTS IN CIFTI FORMAT
##
## DONE!
bglmL <- BayesGLM_cifti(
paste0("data/MSC03_motor07_run0", seq(2), ".4k.dtseries.nii"),
brainstructures="left", surfL = load_surf("left", resamp_res=4000),
onsets=onsets, ar_order=0, EM=FALSE, Bayes=FALSE, resamp_res=NULL, TR=2.2
)
##
## SETTING UP DATA
## .. reading in data for session 1
## .. reading in data for session 2
## MAKING DESIGN MATRICES
## RUNNING MODELS
##
## .. LEFT CORTEX ANALYSIS
## 332 locations removed due to NA or NaN values in at least one scan.
## PUTTING RESULTS IN CIFTI FORMAT
##
## DONE!
combine_and_plot_GLMThis function should combine classical, multi-session GLM results from the left and right hemispheres, and then plot the result. To make it work, add a line that uses combine_xifti to merge the left and right hemisphere "xifti" objects.
combine_and_plot_GLM <- function(
left_GLM, right_GLM,
what=c("estimates", "SE_estimates", "resids"),
...){
# Argument checks.
stopifnot(is.xifti(left_GLM$betas_classical$session_avg))
stopifnot(is.xifti(right_GLM$betas_classical$session_avg))
what <- match.arg(what, c("estimates", "SE_estimates", "resids"))
# Get separate `"xifti"` objects.
if (what=="estimates"){
left_GLM <- left_GLM$betas_classical$session_avg
right_GLM <- right_GLM$betas_classical$session_avg
} else if (what=="SE_estimates") {
# Left
left_mask <- left_GLM$GLMs_classical$cortexL$session_avg$mask
x <- left_GLM$GLMs_classical$cortexL$session_avg$SE_estimates
left_GLM <- newdata_xifti(
left_GLM$betas_classical$session_avg,
x[left_mask,,drop=FALSE]
)
# Right
right_mask <- right_GLM$GLMs_classical$cortexR$session_avg$mask
x <- right_GLM$GLMs_classical$cortexR$session_avg$SE_estimates
right_GLM <- newdata_xifti(
right_GLM$betas_classical$session_avg,
x[right_mask,,drop=FALSE]
)
} else if (what=="resids") {
left_GLM <- newdata_xifti(
left_GLM$betas_classical$session_avg,
left_GLM$GLMs_classical$cortexL$session_avg$resids
)
right_GLM <- newdata_xifti(
right_GLM$betas_classical$session_avg,
right_GLM$GLMs_classical$cortexR$session_avg$resids
)
} else { stop() }
# [Use `combine_xifti` here.]
both_GLM <- combine_xifti(left_GLM, right_GLM)
plot(both_GLM, ...)
}
Try it out!
combine_and_plot_GLM(bglmL, bglmR, "SE")
## `zlim` not provided: using color range 0 - 0.598 (data limits: 0.0725 - 3.58).
Lastly, let’s check out an example of using view_comp. We will plot the GLM result alongside its activations calculated by statistical inference.
# Write and save the GLM result plot.
p1 <- plot(bglmR, fname="glmR_vals", title="GLM result", together="leg")
# Write and save the activations plot.
p2 <- plot(id_activations_cifti(bglmR)$activations_xifti, fname="glmR_act", title="Thresholded activation", together="leg")
# Composite
view_comp(c(p1, p2), nrow=1)
Left foot task activation (right cortex)
Read in cortex CIFTI data without using read_xifti! To work around this challenge, use separate_cifti to split the CIFTI file into GIFTI files, and then read those in with read_xifti2. Check your work by comparing the result to read_xifti. Extra bonus if you can find a way to obtain the subcortical data too.
xii1 <- read_xifti("data/MSC03_motor07_run01.4k.dtseries.nii")
xii2_sep <- separate_cifti("data/MSC03_motor07_run01.4k.dtseries.nii")
xii2 <- read_xifti2(
xii2_sep["cortexL"], xii2_sep["ROIcortexL"],
xii2_sep["cortexR"], xii2_sep["ROIcortexR"]
)
stopifnot(max(as.matrix(xii1-xii2)) == 0)
# Extra bonus
xii1 <- read_xifti(
"data/MSC03_motor07_run01.4k.dtseries.nii",
brainstructures="all"
)
xii2_sep <- separate_cifti(
"data/MSC03_motor07_run01.4k.dtseries.nii",
brainstructures="all"
)
xii2_ctx <- read_xifti2(
xii2_sep["cortexL"], xii2_sep["ROIcortexL"],
xii2_sep["cortexR"], xii2_sep["ROIcortexR"]
)
xii2_sub <- as.xifti(
subcortVol = RNifti::readNifti(xii2_sep["subcortVol"]),
subcortLabs = RNifti::readNifti(xii2_sep["subcortLabs"]),
subcortMask = RNifti::readNifti(xii2_sep["ROIsubcortVol"])
)
xii2 <- combine_xifti(xii2_ctx, xii2_sub)
stopifnot(max(as.matrix(xii1-xii2)) == 0)
ciftiTools can create high-quality plot suitable for publication. Here are a few tips:
fname rather than taking screenshots of the OpenGL window.width and height to increase the plot size. The limits are equal to your screen’s resolution.view_comp can create composite images, but it’s also possible to put together plots by importing them into an image-editing software like GIMP, or to create subplots within LaTeX. For this same reason, you may want to omit title and instead add text above or below the plot later. You can also use legend_embed=FALSE to composite the legend yourself.dlabel files where colors are repeated for different keys, you can use borders=TRUE to delineate them better on cortex plots.material=list(smooth=FALSE, lit=FALSE) to avoid smoothing colors and shading, allowing precise selection of region by color/value for further editing.Please see the Appendix C in the ciftiTools paper for guidance on embedding ciftiTools plots in R Markdown documents.